Disclaimer

This work is very preliminary as I get back into the coding swing of things. Data wrangling and figure generation will be done via R, but the rest of the project will be done using good ol’ microsoft products. This is just an entry point into data crunching and should by no means be considered a final product.

Also, I’m not great at this but whatever. I could automate this, but I’ll figure that out shortly!

Methodology

SNOTEL data was provided by the NRCS. Data was cleaned by removing outliers that are likely implausible; any year with more than 15 observations missing was removed. Temperatures were adjusted using the Morrisey method for stations identified by Ma et al (2019) due to SNOTEL temperature sensor changes, with the adjustment applied to pre-sensor change data. Daily mean observations were detrended to determine whether values were increasing or decreasing from the entire time series trend. Daily mean temperatures were first averaged by water year, with all water year means then averaged by day of water year. The mean temperature by day for the period of record was averaged. To find the standard deviation, the daily mean temperatures by water year was subtracted from the averaged mean temperature by day for the period of record. All water year means averaged by day of water year were subtracted from the temperature mean. The resulting values were then added together to find the “residual” of the daily mean temperatures by water year. The standard deviation was then computed from those residuals, with trends analyzed by Mann‐Kendall significance test and Theil‐Sen’s rate of change. Significant trends are identified with p-values of less than 0.10.

Morrisey Method

The Morrisey Method is taken from Ma, Fassnacht and Kampf..

In R script: T(adjusted) = 5.3x10(-7)xT(old)4+3.72x10(-5)xT(old)3-2.16x10(-3)xT(old)2-7.32x10^(-2)xT(old)+1.37

In the Ma et al. spreadsheet, H1 is Morrisey, H2 is Oiler

San Juan Area SNOTEL sites:

Beartown 327 Original 4/18/2005

Cascade #2 389 Morrisey 6/18/2004

Cumbres Trestle 431 Oyler -> Morrisey 6/8/2005

Idarado 538 Morrisey 7/12/2004

Lone Cone 589 Oyler -> Morrisey 6/22/2005

Middle Creek 624 Morrisey 6/26/2006

Mineral Creek 629 Oyler ?? -> Morrisey 6/22/2004

Molas Lake 632 NSCE-Morrisey, Bias-Original 10/2/2003

Red Mountain Pass 713 Morrisey 8/18/2004

Scotch Creek 739 Morrisey 6/15/2004

Slumgullion 762 Morrisey 6/26/2006

Spud Mountain 780 Morrisey 6/28/2004

Stump Lakes 797 Morrisey 7/22/2005

Upper Rio Grande 839 Oyler -> Morrisey 5/26/2004 *Why is this in red?

Upper San Juan 840 Morrisey 4/7/2004

Vallecito 843 Morrisey 7/22/2005

Wolf Creek Summit 874 Morrisey 7/12/2004

Data from thesis research:


SNOTEL_san_juan_Area <- snotel_download(site_id = c(327, 389, 431, 538, 589, 624, 629, 632, 713, 739, 762, 780, 797, 839, 840, 843, 874), path = tempdir('../data'), internal = TRUE)

write.csv(SNOTEL_san_juan_Area,"C:/Users/13074/Documents/ESS580/thesis_project/San_Juan_area/data_raw/snotel_san_juans.csv", row.names = FALSE) #write in the raw data

Beartown 327

Original

SNOTEL_san_juan_Area <- read.csv("C:/Users/13074/Documents/ESS580/thesis_project/San_Juan_area/data_raw/snotel_san_juans.csv", header = TRUE)

snotel_327 <- SNOTEL_san_juan_Area %>% 
  filter(site_id == "327")
#str(snotel_327) # check the date, usually a character.  

snotel_327$Date <- as.Date(snotel_327$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_327_clean <- snotel_327 %>% # filter for the timeframe
  filter(Date >= "1983-08-10" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_327_clean <- snotel_327_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_327_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_327_clean <- snotel_327_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_min > -40)
ggplot(snotel_327_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
snotel_327_cull_count <- snotel_327_clean %>% 
  filter(temperature_min < -40) %>% 
  count(waterYear)

snotel_327_cull_count
## # A tibble: 0 x 2
## # Groups:   waterYear [0]
## # ... with 2 variables: waterYear <dbl>, n <int>
# filtering for too few observations in a year
snotel_327_cull_count_days <- snotel_327_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_327_cull_count_days
## # A tibble: 3 x 2
## # Groups:   waterYear [3]
##   waterYear     n
##       <dbl> <int>
## 1      1983    52
## 2      1994   337
## 3      2003   346

1983, 1994, 2003 need to be culled.

snotel_327_clean_culled <- snotel_327_clean %>% 
  filter(waterYear != "1983" & waterYear != "1994" & waterYear != "2003") #%>% 
  #filter(temperature_mean > -49)

Beartown (327) NOT ADJUSTED.

2005-04-18 installed ext range sensor. Original

snotel_327_adjusted <- snotel_327_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2005-04-18", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

327 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_327 <- snotel_327_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_327 <- yearly_wy_aver_327 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_327 <- daily_wy_aver_327 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_327$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_327 <-daily_wy_aver_327 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_327$date_temp <- signif(daily_wy_aver2_327$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_327, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

327 SD

standard_dev_327 <- daily_wy_aver_327 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_327 <- standard_dev_327 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_327 <- standard_dev_all_327 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_327 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 8.215357
1985 8.122497
1986 7.129599
1987 7.619777
1988 8.257318
1989 8.044976
1990 7.675963
1991 7.684970
1992 7.079517
1993 7.642300
1995 7.304999
1996 7.681058
1997 7.711950
1998 8.105700
1999 6.825322
2000 7.526413
2001 7.912097
2002 8.099855
2004 7.622635
2005 7.225545
2006 7.107924
2007 7.426444
2008 7.802192
2009 6.893765
2010 7.820347
2011 7.632307
2012 7.466342
2013 8.000100
2014 7.260505
2015 6.671562
2016 7.474460
2017 7.072348
2018 6.933787
2019 7.658054
2020 7.769482
2021 7.841403
2022 7.318492
ggplot(standard_dev_all_327, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 327 average temperatures for water years 2005-2021

MK & SS for 327 (non-corrected)

sd_mk_327 <- mk.test(standard_dev_all_327$sd_2)
print(sd_mk_327)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_327$sd_2
## z = -1.818, n = 37, p-value = 0.06907
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## -140.0000000 5846.0000000   -0.2102102
sd_sens_327 <- sens.slope(standard_dev_all_327$sd_2)
print(sd_sens_327)
## 
##  Sen's slope
## 
## data:  standard_dev_all_327$sd_2
## z = -1.818, n = 37, p-value = 0.06907
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.025461316  0.001275883
## sample estimates:
## Sen's slope 
## -0.01223021

Cascade #2 389

Morrisey 6/18/2004

snotel_389 <- SNOTEL_san_juan_Area %>% 
  filter(site_id == "389")
#str(snotel_389) # check the date, usually a character.  

snotel_389$Date <- as.Date(snotel_389$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_389_clean <- snotel_389 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_389_clean <- snotel_389_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_389_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_389_clean <- snotel_389_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_min > -40) %>% 
  filter(temperature_mean > -40)
ggplot(snotel_389_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

temp_389_xts <- xts(snotel_389_clean$temperature_mean, order.by = snotel_389_clean$Date)

dygraph(temp_389_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
snotel_389_cull_count <- snotel_389_clean %>% 
  filter(temperature_min > -40) %>% 
  count(waterYear)

snotel_389_cull_count
## # A tibble: 40 x 2
## # Groups:   waterYear [40]
##    waterYear     n
##        <dbl> <int>
##  1      1983   325
##  2      1984   350
##  3      1985   350
##  4      1986   365
##  5      1987   358
##  6      1988   366
##  7      1989   330
##  8      1990   348
##  9      1991   364
## 10      1992   366
## # ... with 30 more rows
# filtering for too few observations in a year
snotel_389_cull_count_days <- snotel_389_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_389_cull_count_days
## # A tibble: 9 x 2
## # Groups:   waterYear [9]
##   waterYear     n
##       <dbl> <int>
## 1      1983   325
## 2      1989   330
## 3      1990   348
## 4      1994   310
## 5      2006   330
## 6      2010   334
## 7      2014   314
## 8      2015   185
## 9      2016   334

1983, 1989, 1990, 1994, 2006, 2010, 2014, 2015, 2016 need to be culled.

snotel_389_clean_culled <- snotel_389_clean %>% 
  filter(waterYear != "1983" & waterYear != "1990" & waterYear != "1994" & waterYear != "2006" & waterYear != "2010" & waterYear != "2014" & waterYear != "2015" & waterYear != "2016") #%>% 
  #filter(temperature_mean > -49)

Cascade #2 (389) NOT ADJUSTED.

Morrisey 6/18/2004

snotel_389_adjusted <- snotel_389_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2004-06-18", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

389 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_389 <- snotel_389_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_389 <- yearly_wy_aver_389 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_389 <- daily_wy_aver_389 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_389$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_389 <-daily_wy_aver_389 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_389$date_temp <- signif(daily_wy_aver2_389$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_389, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

389 SD

standard_dev_389 <- daily_wy_aver_389 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_389 <- standard_dev_389 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_389 <- standard_dev_all_389 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_389 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 10.776321
1985 7.879611
1986 10.887499
1987 9.408888
1988 10.711609
1989 10.418326
1991 10.101153
1992 8.637142
1993 8.884440
1995 9.371322
1996 9.429341
1997 10.027003
1998 9.698491
1999 9.029327
2000 9.288091
2001 10.630687
2002 10.346726
2003 9.727202
2004 9.282545
2005 8.472227
2007 9.686674
2008 9.541070
2009 8.659960
2011 9.334867
2012 9.418612
2013 9.751020
2017 9.144382
2018 8.911243
2019 9.427120
2020 9.880776
2021 9.863539
2022 9.441803
ggplot(standard_dev_all_389, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 389 average temperatures for water years 2005-2021

MK & SS for 389 (non-corrected)

sd_mk_389 <- mk.test(standard_dev_all_389$sd_2)
print(sd_mk_389)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_389$sd_2
## z = -0.95677, n = 32, p-value = 0.3387
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  -60.0000000 3802.6666667   -0.1209677
sd_sens_389 <- sens.slope(standard_dev_all_389$sd_2)
print(sd_sens_389)
## 
##  Sen's slope
## 
## data:  standard_dev_all_389$sd_2
## z = -0.95677, n = 32, p-value = 0.3387
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.04702986  0.01336703
## sample estimates:
## Sen's slope 
## -0.01646233

Corrected

389 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_389_ad <- snotel_389_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_389_ad <- yearly_wy_aver_389_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_389_ad <- daily_wy_aver_389_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_389_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_389_ad <-daily_wy_aver_389_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_389_ad$date_temp_ad <- signif(daily_wy_aver2_389_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_389_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

389 SS (corrected)

standard_dev_389_ad <- daily_wy_aver_389_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_389_ad <- standard_dev_389_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_389_ad <- standard_dev_all_389_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_389_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 10.061639
1985 7.175032
1986 9.930924
1987 8.682951
1988 9.868631
1989 9.733860
1991 9.472969
1992 8.002025
1993 8.283300
1995 8.657379
1996 8.697682
1997 9.311693
1998 8.929771
1999 8.328247
2000 8.510889
2001 9.808164
2002 9.563945
2003 8.940614
2004 8.587303
2005 8.472294
2007 9.686174
2008 9.541322
2009 8.659292
2011 9.334237
2012 9.418205
2013 9.750812
2017 9.144257
2018 8.911176
2019 9.426664
2020 9.880403
2021 9.863015
2022 9.441185
ggplot(standard_dev_all_389_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 389 average temperatures for water years 1986-2021

MK & SS 389 (corrected)

sd_mk_389_ad <- mk.test(standard_dev_all_389_ad$sd_2)
print(sd_mk_389_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_389_ad$sd_2
## z = 1.0216, n = 32, p-value = 0.307
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##   64.0000000 3802.6666667    0.1290323
sd_sens_389_ad <- sens.slope(standard_dev_all_389_ad$sd_2)
print(sd_sens_389_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_389_ad$sd_2
## z = 1.0216, n = 32, p-value = 0.307
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01140359  0.04896578
## sample estimates:
## Sen's slope 
##  0.01785583

Cumbres Trestle 431

Oyler -> Morrisey 6/8/2005

snotel_431 <- SNOTEL_san_juan_Area %>% 
  filter(site_id == "431")
#str(snotel_431) # check the date, usually a character.  

snotel_431$Date <- as.Date(snotel_431$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_431_clean <- snotel_431 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_431_clean <- snotel_431_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_431_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_431_clean <- snotel_431_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_min > -40) %>% 
  filter(temperature_min < 25)
ggplot(snotel_431_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
snotel_431_cull_count <- snotel_431_clean %>% 
  filter(temperature_min > -40) %>% 
  count(waterYear)

snotel_431_cull_count
## # A tibble: 36 x 2
## # Groups:   waterYear [36]
##    waterYear     n
##        <dbl> <int>
##  1      1987     1
##  2      1988   341
##  3      1989   364
##  4      1990   348
##  5      1991   341
##  6      1992   312
##  7      1993   356
##  8      1994   317
##  9      1995   364
## 10      1996   364
## # ... with 26 more rows
# filtering for too few observations in a year
snotel_431_cull_count_days <- snotel_431_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_431_cull_count_days
## # A tibble: 7 x 2
## # Groups:   waterYear [7]
##   waterYear     n
##       <dbl> <int>
## 1      1987     1
## 2      1988   341
## 3      1990   348
## 4      1991   341
## 5      1992   312
## 6      1994   317
## 7      2014   288

1987, 1988, 1990, 1991, 1992, 1994, 2014 need to be culled.

snotel_431_clean_culled <- snotel_431_clean %>% 
  filter(waterYear != "1987" & waterYear != "1988" & waterYear != "1990" & waterYear != "1991" & waterYear != "1992" & waterYear != "1994" & waterYear != "2014") #%>% 
  #filter(temperature_mean > -49)
ggplot(snotel_431_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_431_xts <- xts(snotel_431_clean_culled$temperature_mean, order.by = snotel_431_clean_culled$Date)

dygraph(temp_431_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
snotel_431_clean_culled <- snotel_431_clean_culled %>% 
  filter(temperature_mean < 20)

temp_431_xts <- xts(snotel_431_clean_culled$temperature_mean, order.by = snotel_431_clean_culled$Date)

dygraph(temp_431_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 

Cumbres Trestle 431

Oyler -> Morrisey 6/8/2005

snotel_431_adjusted <- snotel_431_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2005-06-08", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

431 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_431 <- snotel_431_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_431 <- yearly_wy_aver_431 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_431 <- daily_wy_aver_431 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_431$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_431 <-daily_wy_aver_431 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_431$date_temp <- signif(daily_wy_aver2_431$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_431, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

431 SD

standard_dev_431 <- daily_wy_aver_431 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_431 <- standard_dev_431 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_431 <- standard_dev_all_431 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_431 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1989 8.484011
1993 8.101376
1995 7.750446
1996 7.952021
1997 8.217323
1998 8.600343
1999 7.226028
2000 8.164156
2001 8.669911
2002 9.005161
2003 8.459072
2004 8.306872
2005 8.304679
2006 7.417238
2007 7.937180
2008 8.259566
2009 7.296965
2010 8.418225
2011 8.231369
2012 8.081428
2013 8.537331
2015 7.219164
2016 7.791453
2017 7.322790
2018 7.113402
2019 7.880447
2020 8.046531
2021 8.095523
2022 7.621146
ggplot(standard_dev_all_431, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 431 average temperatures for water years 2005-2021

MK & SS for 431 (non-corrected)

sd_mk_431 <- mk.test(standard_dev_all_431$sd_2)
print(sd_mk_431)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_431$sd_2
## z = -1.9321, n = 29, p-value = 0.05335
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## -104.0000000 2842.0000000   -0.2561576
sd_sens_431 <- sens.slope(standard_dev_all_431$sd_2)
print(sd_sens_431)
## 
##  Sen's slope
## 
## data:  standard_dev_all_431$sd_2
## z = -1.9321, n = 29, p-value = 0.05335
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.041299369  0.001003281
## sample estimates:
## Sen's slope 
## -0.01890302

Corrected

431 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_431_ad <- snotel_431_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_431_ad <- yearly_wy_aver_431_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_431_ad <- daily_wy_aver_431_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_431_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_431_ad <-daily_wy_aver_431_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_431_ad$date_temp_ad <- signif(daily_wy_aver2_431_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_431_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

431 SS (corrected)

standard_dev_431_ad <- daily_wy_aver_431_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_431_ad <- standard_dev_431_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_431_ad <- standard_dev_all_431_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_431_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1989 7.914872
1993 7.518415
1995 7.173944
1996 7.365597
1997 7.634730
1998 7.970563
1999 6.682639
2000 7.537564
2001 8.039524
2002 8.344015
2003 7.813684
2004 7.737034
2005 7.599212
2006 7.418573
2007 7.938448
2008 8.260644
2009 7.298209
2010 8.419846
2011 8.232506
2012 8.082571
2013 8.538590
2015 7.220737
2016 7.793079
2017 7.324026
2018 7.114853
2019 7.881605
2020 8.047598
2021 8.096826
2022 7.622159
ggplot(standard_dev_all_431_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 431 average temperatures for water years 1986-2021

MK & SS 431 (corrected)

sd_mk_431_ad <- mk.test(standard_dev_all_431_ad$sd_2)
print(sd_mk_431_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_431_ad$sd_2
## z = 0.99418, n = 29, p-value = 0.3201
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##   54.0000000 2842.0000000    0.1330049
sd_sens_431_ad <- sens.slope(standard_dev_all_431_ad$sd_2)
print(sd_sens_431_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_431_ad$sd_2
## z = 0.99418, n = 29, p-value = 0.3201
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01334615  0.03134200
## sample estimates:
## Sen's slope 
##  0.01029449

Idarado 538

Morrisey 7/12/2004

snotel_538 <- SNOTEL_san_juan_Area %>% 
  filter(site_id == "538")
#str(snotel_538) # check the date, usually a character.  

snotel_538$Date <- as.Date(snotel_538$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_538_clean <- snotel_538 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_538_clean <- snotel_538_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_538_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_538_clean <- snotel_538_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_mean > -40)# %>% 
  #filter(temperature_min < 25)
ggplot(snotel_538_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
snotel_538_cull_count <- snotel_538_clean %>% 
  filter(temperature_min > -40) %>% 
  count(waterYear)

snotel_538_cull_count
## # A tibble: 36 x 2
## # Groups:   waterYear [36]
##    waterYear     n
##        <dbl> <int>
##  1      1987   355
##  2      1988   366
##  3      1989   364
##  4      1990   365
##  5      1991   337
##  6      1992   358
##  7      1993   346
##  8      1994   351
##  9      1995   364
## 10      1996   365
## # ... with 26 more rows
# filtering for too few observations in a year
snotel_538_cull_count_days <- snotel_538_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_538_cull_count_days
## # A tibble: 2 x 2
## # Groups:   waterYear [2]
##   waterYear     n
##       <dbl> <int>
## 1      1991   337
## 2      1993   346

1991 and 1993 need to be culled.

snotel_538_clean_culled <- snotel_538_clean %>% 
  filter(waterYear != "1991" & waterYear != "1993") # & waterYear != "1990" & waterYear != "1991" & waterYear != "1992" & waterYear != "1994" & waterYear != "2014") #%>% 
  #filter(temperature_mean > -49)
ggplot(snotel_538_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_538_xts <- xts(snotel_538_clean_culled$temperature_mean, order.by = snotel_538_clean_culled$Date)

dygraph(temp_538_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
snotel_538_clean_culled <- snotel_538_clean_culled %>% 
  filter(temperature_mean < 20)

temp_538_xts <- xts(snotel_538_clean_culled$temperature_mean, order.by = snotel_538_clean_culled$Date)

dygraph(temp_538_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 

Idarado 538

Morrisey 7/12/2004

snotel_538_adjusted <- snotel_538_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2004-07-12", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

538 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_538 <- snotel_538_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_538 <- yearly_wy_aver_538 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_538 <- daily_wy_aver_538 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_538$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_538 <-daily_wy_aver_538 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_538$date_temp <- signif(daily_wy_aver2_538$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_538, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

538 SD

standard_dev_538 <- daily_wy_aver_538 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_538 <- standard_dev_538 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_538 <- standard_dev_all_538 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_538 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 8.006364
1988 8.601864
1989 8.632529
1990 8.392986
1992 7.829013
1994 10.169815
1995 7.658097
1996 7.925011
1997 7.934668
1998 8.455869
1999 7.252626
2000 7.992580
2001 8.428201
2002 8.883209
2003 8.309823
2004 8.066471
2005 7.228205
2006 7.518896
2007 7.890984
2008 8.117130
2009 7.048527
2010 8.062733
2011 7.859864
2012 7.787906
2013 8.410412
2014 7.465408
2015 6.906228
2016 7.611126
2017 7.235848
2018 7.239780
2019 7.718304
2020 7.976082
2021 8.027478
2022 7.664387
ggplot(standard_dev_all_538, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 538 average temperatures for water years 2005-2021

MK & SS for 538 (non-corrected)

sd_mk_538 <- mk.test(standard_dev_all_538$sd_2)
print(sd_mk_538)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_538$sd_2
## z = -2.6684, n = 34, p-value = 0.007621
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## -181.0000000 4550.3333333   -0.3226381
sd_sens_538 <- sens.slope(standard_dev_all_538$sd_2)
print(sd_sens_538)
## 
##  Sen's slope
## 
## data:  standard_dev_all_538$sd_2
## z = -2.6684, n = 34, p-value = 0.007621
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.042981556 -0.006896251
## sample estimates:
## Sen's slope 
## -0.02384088

Corrected

538 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_538_ad <- snotel_538_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_538_ad <- yearly_wy_aver_538_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_538_ad <- daily_wy_aver_538_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_538_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_538_ad <-daily_wy_aver_538_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_538_ad$date_temp_ad <- signif(daily_wy_aver2_538_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_538_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

538 SS (corrected)

standard_dev_538_ad <- daily_wy_aver_538_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_538_ad <- standard_dev_538_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_538_ad <- standard_dev_all_538_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_538_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 7.447584
1988 8.019241
1989 8.069873
1990 7.812557
1992 7.272544
1994 9.549437
1995 7.113641
1996 7.353811
1997 7.396458
1998 7.860736
1999 6.739299
2000 7.408528
2001 7.835100
2002 8.260161
2003 7.706375
2004 7.450714
2005 7.228187
2006 7.519259
2007 7.891254
2008 8.117349
2009 7.048722
2010 8.062844
2011 7.860215
2012 7.788573
2013 8.411054
2014 7.465688
2015 6.906947
2016 7.611742
2017 7.235817
2018 7.240024
2019 7.718445
2020 7.976612
2021 8.027693
2022 7.664991
ggplot(standard_dev_all_538_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 538 average temperatures for water years 1986-2021

MK & SS 538 (corrected)

sd_mk_538_ad <- mk.test(standard_dev_all_538_ad$sd_2)
print(sd_mk_538_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_538_ad$sd_2
## z = 0.14824, n = 34, p-value = 0.8821
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 1.100000e+01 4.550333e+03 1.960784e-02
sd_sens_538_ad <- sens.slope(standard_dev_all_538_ad$sd_2)
print(sd_sens_538_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_538_ad$sd_2
## z = 0.14824, n = 34, p-value = 0.8821
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01508416  0.01703536
## sample estimates:
##  Sen's slope 
## 0.0009105465

Lone Cone 589

Oyler -> Morrisey 6/22/2005

snotel_589 <- SNOTEL_san_juan_Area %>% 
  filter(site_id == "589")
#str(snotel_589) # check the date, usually a character.  

snotel_589$Date <- as.Date(snotel_589$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_589_clean <- snotel_589 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_589_clean <- snotel_589_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_589_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_589_clean <- snotel_589_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_mean > -40)# %>% 
  #filter(temperature_min < 25)
ggplot(snotel_589_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
snotel_589_cull_count <- snotel_589_clean %>% 
  filter(temperature_min > -40) %>% 
  count(waterYear)

snotel_589_cull_count
## # A tibble: 37 x 2
## # Groups:   waterYear [37]
##    waterYear     n
##        <dbl> <int>
##  1      1986     1
##  2      1987   358
##  3      1988   366
##  4      1989   362
##  5      1990   365
##  6      1991   359
##  7      1992   366
##  8      1993   350
##  9      1994   337
## 10      1995   364
## # ... with 27 more rows
# filtering for too few observations in a year
snotel_589_cull_count_days <- snotel_589_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_589_cull_count_days
## # A tibble: 10 x 2
## # Groups:   waterYear [10]
##    waterYear     n
##        <dbl> <int>
##  1      1986     1
##  2      1994   337
##  3      1997   346
##  4      2006   317
##  5      2009   314
##  6      2010   332
##  7      2011   309
##  8      2012   322
##  9      2013   331
## 10      2014   335

1986, 1994, 1997, 2006, 2009, 2010, 2011, 2012, 2013, 2014 need to be culled.

snotel_589_clean_culled <- snotel_589_clean %>% 
  filter(waterYear != "1986" & waterYear != "1994" & waterYear != "1997" & waterYear != "2006" & waterYear != "2009" & waterYear != "2010" & waterYear != "2011" & waterYear != "2012" & waterYear != "2013" & waterYear != "2014") #%>% 
  #filter(temperature_mean > -49)
ggplot(snotel_589_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_589_xts <- xts(snotel_589_clean_culled$temperature_mean, order.by = snotel_589_clean_culled$Date)

dygraph(temp_589_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
snotel_589_clean_culled <- snotel_589_clean_culled %>% 
  filter(temperature_mean < 20)

temp_589_xts <- xts(snotel_589_clean_culled$temperature_mean, order.by = snotel_589_clean_culled$Date)

dygraph(temp_589_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 

Lone Cone 589

Oyler -> Morrisey 6/22/2005

snotel_589_adjusted <- snotel_589_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2005-06-22", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

589 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_589 <- snotel_589_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_589 <- yearly_wy_aver_589 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_589 <- daily_wy_aver_589 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_589$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_589 <-daily_wy_aver_589 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_589$date_temp <- signif(daily_wy_aver2_589$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_589, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

589 SD

standard_dev_589 <- daily_wy_aver_589 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_589 <- standard_dev_589 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_589 <- standard_dev_all_589 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_589 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 8.224944
1988 8.773281
1989 8.776693
1990 8.580289
1991 8.502502
1992 7.975611
1993 8.199081
1995 7.946318
1996 8.231469
1998 8.706109
1999 7.431979
2000 8.232988
2001 8.771085
2002 9.082003
2003 8.732769
2004 8.380419
2005 8.458900
2007 8.212611
2008 8.408220
2015 7.103666
2016 7.842413
2017 7.527960
2018 7.506884
2019 8.029199
2020 8.271099
2021 8.278625
2022 8.020852
ggplot(standard_dev_all_589, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 589 average temperatures for water years 2005-2021

MK & SS for 589 (non-corrected)

sd_mk_589 <- mk.test(standard_dev_all_589$sd_2)
print(sd_mk_589)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_589$sd_2
## z = -1.8345, n = 27, p-value = 0.06658
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  -89.0000000 2301.0000000   -0.2535613
sd_sens_589 <- sens.slope(standard_dev_all_589$sd_2)
print(sd_sens_589)
## 
##  Sen's slope
## 
## data:  standard_dev_all_589$sd_2
## z = -1.8345, n = 27, p-value = 0.06658
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.043902404  0.001230005
## sample estimates:
## Sen's slope 
## -0.02183399

Corrected

589 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_589_ad <- snotel_589_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_589_ad <- yearly_wy_aver_589_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_589_ad <- daily_wy_aver_589_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_589_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_589_ad <-daily_wy_aver_589_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_589_ad$date_temp_ad <- signif(daily_wy_aver2_589_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_589_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

589 SS (corrected)

standard_dev_589_ad <- daily_wy_aver_589_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_589_ad <- standard_dev_589_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_589_ad <- standard_dev_all_589_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_589_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 7.671688
1988 8.200568
1989 8.222864
1990 7.999469
1991 7.962227
1992 7.420415
1993 7.639311
1995 7.378994
1996 7.641722
1998 8.086353
1999 6.914950
2000 7.633161
2001 8.156564
2002 8.440080
2003 8.103839
2004 7.827141
2005 7.781158
2007 8.212720
2008 8.408841
2015 7.104293
2016 7.842763
2017 7.527954
2018 7.507329
2019 8.029667
2020 8.271310
2021 8.278726
2022 8.021229
ggplot(standard_dev_all_589_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 589 average temperatures for water years 1986-2021

MK & SS 589 (corrected)

sd_mk_589_ad <- mk.test(standard_dev_all_589_ad$sd_2)
print(sd_mk_589_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_589_ad$sd_2
## z = 0.6671, n = 27, p-value = 0.5047
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 3.300000e+01 2.301000e+03 9.401709e-02
sd_sens_589_ad <- sens.slope(standard_dev_all_589_ad$sd_2)
print(sd_sens_589_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_589_ad$sd_2
## z = 0.6671, n = 27, p-value = 0.5047
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01153608  0.02498422
## sample estimates:
## Sen's slope 
## 0.008250776

Middle Creek 624

Morrisey 6/26/2006

snotel_624 <- SNOTEL_san_juan_Area %>% 
  filter(site_id == "624")
#str(snotel_624) # check the date, usually a character.  

snotel_624$Date <- as.Date(snotel_624$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_624_clean <- snotel_624 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_624_clean <- snotel_624_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_624_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_624_clean <- snotel_624_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  filter(temperature_mean > -40) %>% 
  filter(temperature_mean < 25)
ggplot(snotel_624_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
snotel_624_cull_count <- snotel_624_clean %>% 
  filter(temperature_min > -40) %>% 
  count(waterYear)

snotel_624_cull_count
## # A tibble: 43 x 2
## # Groups:   waterYear [43]
##    waterYear     n
##        <dbl> <int>
##  1      1980    59
##  2      1981   354
##  3      1982     1
##  4      1983   312
##  5      1984   360
##  6      1985   362
##  7      1986   364
##  8      1987   360
##  9      1988   340
## 10      1989   363
## # ... with 33 more rows
# filtering for too few observations in a year
snotel_624_cull_count_days <- snotel_624_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_624_cull_count_days
## # A tibble: 9 x 2
## # Groups:   waterYear [9]
##   waterYear     n
##       <dbl> <int>
## 1      1980    59
## 2      1982     1
## 3      1983   312
## 4      1988   340
## 5      1991   338
## 6      1993   330
## 7      1994   330
## 8      2021   338
## 9      2022   348

1980, 1982, 1983, 1988, 1991, 1993, 1994, 2021, 2022 need to be culled.

snotel_624_clean_culled <- snotel_624_clean %>% 
  filter(waterYear != "1980" & waterYear != "1981" & waterYear != "1982" & waterYear != "1983" & waterYear != "1988" & waterYear != "1991" & waterYear != "1993" & waterYear != "1994" & waterYear != "2021" & waterYear != "2022")# & waterYear != "2014") #%>% 
  #filter(temperature_mean > -49)

ggplot(snotel_624_clean_culled, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

ggplot(snotel_624_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_624_xts <- xts(snotel_624_clean_culled$temperature_mean, order.by = snotel_624_clean_culled$Date)

dygraph(temp_624_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
snotel_624_clean_culled <- snotel_624_clean_culled %>% 
  filter(temperature_mean < 20)

temp_624_xts <- xts(snotel_624_clean_culled$temperature_mean, order.by = snotel_624_clean_culled$Date)

dygraph(temp_624_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 

Middle Creek 624 Morrisey 6/26/2006

snotel_624_adjusted <- snotel_624_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2006-06-26", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

624 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_624 <- snotel_624_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_624 <- yearly_wy_aver_624 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_624 <- daily_wy_aver_624 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_624$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_624 <-daily_wy_aver_624 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_624$date_temp <- signif(daily_wy_aver2_624$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_624, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

624 SD

standard_dev_624 <- daily_wy_aver_624 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_624 <- standard_dev_624 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_624 <- standard_dev_all_624 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_624 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 8.239005
1985 8.207104
1986 7.191318
1987 7.871848
1989 8.454185
1990 8.200220
1992 6.438238
1995 7.565327
1996 7.835558
1997 7.996219
1998 8.415139
1999 7.174609
2000 7.839800
2001 8.355597
2002 8.492066
2003 8.150297
2004 7.912597
2005 7.915951
2006 7.936724
2007 7.640492
2008 8.021623
2009 7.171130
2010 8.259271
2011 7.932618
2012 7.791063
2013 8.259320
2014 7.372950
2015 6.806548
2016 7.548199
2017 7.108731
2018 7.036935
2019 7.675916
2020 7.741703
ggplot(standard_dev_all_624, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 624 average temperatures for water years 2005-2021

MK & SS for 624 (non-corrected)

sd_mk_624 <- mk.test(standard_dev_all_624$sd_2)
print(sd_mk_624)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_624$sd_2
## z = -1.9678, n = 33, p-value = 0.04909
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## -128.0000000 4165.3333333   -0.2424242
sd_sens_624 <- sens.slope(standard_dev_all_624$sd_2)
print(sd_sens_624)
## 
##  Sen's slope
## 
## data:  standard_dev_all_624$sd_2
## z = -1.9678, n = 33, p-value = 0.04909
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.0333661824 -0.0003479011
## sample estimates:
## Sen's slope 
## -0.01662145

Corrected

624 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_624_ad <- snotel_624_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_624_ad <- yearly_wy_aver_624_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_624_ad <- daily_wy_aver_624_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_624_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_624_ad <-daily_wy_aver_624_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_624_ad$date_temp_ad <- signif(daily_wy_aver2_624_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_624_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

624 SS (corrected)

standard_dev_624_ad <- daily_wy_aver_624_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_624_ad <- standard_dev_624_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_624_ad <- standard_dev_all_624_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_624_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 7.738801
1985 7.713059
1986 6.725366
1987 7.374899
1989 7.942209
1990 7.687289
1992 6.080980
1995 7.080976
1996 7.336792
1997 7.519955
1998 7.885653
1999 6.718142
2000 7.329553
2001 7.833709
2002 7.944757
2003 7.617514
2004 7.442397
2005 7.390605
2006 7.304030
2007 7.639661
2008 8.021493
2009 7.170919
2010 8.259153
2011 7.932368
2012 7.790490
2013 8.259010
2014 7.372461
2015 6.806595
2016 7.547948
2017 7.107845
2018 7.036548
2019 7.675791
2020 7.741347
ggplot(standard_dev_all_624_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 624 average temperatures for water years 1986-2021

MK & SS 624 (corrected)

sd_mk_624_ad <- mk.test(standard_dev_all_624_ad$sd_2)
print(sd_mk_624_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_624_ad$sd_2
## z = 0.35637, n = 33, p-value = 0.7216
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 2.400000e+01 4.165333e+03 4.545455e-02
sd_sens_624_ad <- sens.slope(standard_dev_all_624_ad$sd_2)
print(sd_sens_624_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_624_ad$sd_2
## z = 0.35637, n = 33, p-value = 0.7216
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01304163  0.02223675
## sample estimates:
## Sen's slope 
## 0.003307866

Molas Lake 632

NSCE-Morrisey, Bias-Original 10/2/2003

snotel_632 <- SNOTEL_san_juan_Area %>% 
  filter(site_id == "632")
#str(snotel_632) # check the date, usually a character.  

snotel_632$Date <- as.Date(snotel_632$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.

#THIS WILL CHANGE FOR EACH STATION
snotel_632_clean <- snotel_632 %>% # filter for the timeframe
  filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
  #filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers   
  addWaterYear() %>% 
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  na.omit()

#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))

snotel_632_clean <- snotel_632_clean %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))
# Check for outliers

ggplot(snotel_632_clean, aes(x = Date, y = temperature_mean)) +
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

snotel_632_clean <- snotel_632_clean %>% 
  mutate(temp_diff = abs(temperature_min - temperature_max)) %>% 
  #filter(temperature_mean > -40) %>% 
  filter(temperature_mean < 25)
ggplot(snotel_632_clean, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”

# filtering for temperature anomalies
snotel_632_cull_count <- snotel_632_clean %>% 
  filter(temperature_min > -40) %>% 
  count(waterYear)

snotel_632_cull_count
## # A tibble: 37 x 2
## # Groups:   waterYear [37]
##    waterYear     n
##        <dbl> <int>
##  1      1986    54
##  2      1987   355
##  3      1988   362
##  4      1989   364
##  5      1990   365
##  6      1991   365
##  7      1992   366
##  8      1993   303
##  9      1994   294
## 10      1995   361
## # ... with 27 more rows
# filtering for too few observations in a year
snotel_632_cull_count_days <- snotel_632_clean %>% 
  group_by(waterYear) %>% 
  count(waterYear) %>% 
  filter(n < 350)

snotel_632_cull_count_days
## # A tibble: 6 x 2
## # Groups:   waterYear [6]
##   waterYear     n
##       <dbl> <int>
## 1      1986    54
## 2      1993   303
## 3      1994   295
## 4      2005   344
## 5      2006   304
## 6      2013   248

1986, 1993, 1994, 2005, 2006, 2013 need to be culled.

snotel_632_clean_culled <- snotel_632_clean %>% 
  filter(waterYear != "1986" & waterYear != "1993" & waterYear != "1994" & waterYear != "2005" & waterYear != "2006" & waterYear != "2013")# & waterYear != "1993" & waterYear != "1994" & waterYear != "2021" & waterYear != "2022") & waterYear != "2014") #%>% 
  #filter(temperature_mean > -49)

ggplot(snotel_632_clean_culled, aes(x = Date, y = temp_diff)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature varience (°C)') + 
  xlab('Date')

ggplot(snotel_632_clean_culled, aes(x = Date, y = temperature_mean)) + 
  geom_point() + #lwd = 2) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

temp_632_xts <- xts(snotel_632_clean_culled$temperature_mean, order.by = snotel_632_clean_culled$Date)

dygraph(temp_632_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 
snotel_632_clean_culled <- snotel_632_clean_culled %>% 
  filter(temperature_mean < 20)

temp_632_xts <- xts(snotel_632_clean_culled$temperature_mean, order.by = snotel_632_clean_culled$Date)

dygraph(temp_632_xts) %>%
  dyAxis("y", label = "Daily mean temperature (°C)") 

Molas Lake 632

NSCE-Morrisey, Bias-Original 10/2/2003

snotel_632_adjusted <- snotel_632_clean_culled %>%
  mutate(temp_ad = if_else(Date < "2003-10-02", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))

Non-corrected

632 Detrended

#using the clean culled df:

#average water year temperature

yearly_wy_aver_632 <- snotel_632_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp = mean(temperature_mean))

#Average temperature by day for all water years:

daily_wy_aver_632 <- yearly_wy_aver_632 %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

daily_wy_aver_632 <- daily_wy_aver_632 %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(daily_wy_aver_632$aver_day_temp))

# try to show all years as means. 
daily_wy_aver2_632 <-daily_wy_aver_632 %>% 
  group_by(waterDay) %>%
  mutate(date_temp = mean(temperature_mean))
  

daily_wy_aver2_632$date_temp <- signif(daily_wy_aver2_632$date_temp,3) #reduce the sig figs


ggplot(daily_wy_aver2_632, aes(x = waterDay, y = date_temp))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

632 SD

standard_dev_632 <- daily_wy_aver_632 %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

standard_dev_all_632 <- standard_dev_632 %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

standard_dev_all_632 <- standard_dev_all_632 %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

standard_dev_all_632 %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 7.613234
1988 8.191410
1989 8.166070
1990 7.646926
1991 7.969691
1992 7.156808
1995 7.489888
1996 7.655841
1997 7.698012
1998 8.157837
1999 6.885991
2000 7.657554
2001 6.788385
2002 8.254532
2003 7.907020
2004 7.308264
2007 7.607147
2008 7.964946
2009 6.957483
2010 7.983696
2011 7.701623
2012 7.490054
2014 7.236529
2015 6.581104
2016 7.409875
2017 7.001441
2018 6.814621
2019 7.582515
2020 7.766470
2021 7.792571
2022 7.357711
ggplot(standard_dev_all_632, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Non-corrected standard deviation of SNOTEL 632 average temperatures for water years 2005-2021

MK & SS for 632 (non-corrected)

sd_mk_632 <- mk.test(standard_dev_all_632$sd_2)
print(sd_mk_632)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_632$sd_2
## z = -1.6317, n = 31, p-value = 0.1028
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  -97.0000000 3461.6666667   -0.2086022
sd_sens_632 <- sens.slope(standard_dev_all_632$sd_2)
print(sd_sens_632)
## 
##  Sen's slope
## 
## data:  standard_dev_all_632$sd_2
## z = -1.6317, n = 31, p-value = 0.1028
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.034372524  0.003422909
## sample estimates:
## Sen's slope 
## -0.01583493

Corrected

632 Detrended (corrected)

#using the clean culled df:
#average water year temperature
yearly_wy_aver_632_ad <- snotel_632_adjusted %>% 
  group_by(waterYear) %>% 
  mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_632_ad <- yearly_wy_aver_632_ad %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_632_ad <- daily_wy_aver_632_ad %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp_ad = mean(daily_wy_aver_632_ad$aver_day_temp_ad))
# try to show all years as means. 
daily_wy_aver2_632_ad <-daily_wy_aver_632_ad %>% 
  group_by(waterDay) %>%
  mutate(date_temp_ad = mean(temp_ad))
  
daily_wy_aver2_632_ad$date_temp_ad <- signif(daily_wy_aver2_632_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_632_ad, aes(x = waterDay, y = date_temp_ad))+
  geom_line(size= 0.7) +
  theme_few() +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

632 SS (corrected)

standard_dev_632_ad <- daily_wy_aver_632_ad %>% 
  group_by(waterYear) %>% 
  #filter(waterYear >= 1987 & waterYear <= 2021) %>% 
  mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>% 
  mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_632_ad <- standard_dev_632_ad %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())
standard_dev_all_632_ad <- standard_dev_all_632_ad %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)
standard_dev_all_632_ad %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1987 7.106416
1988 7.678134
1989 7.662133
1990 7.141729
1991 7.492618
1992 6.671512
1995 6.957618
1996 7.110780
1997 7.178235
1998 7.586041
1999 6.397699
2000 7.103150
2001 6.358660
2002 7.672831
2003 7.338249
2004 7.309361
2007 7.606845
2008 7.965266
2009 6.957636
2010 7.984148
2011 7.701741
2012 7.490245
2014 7.236308
2015 6.581364
2016 7.409996
2017 7.000722
2018 6.815246
2019 7.583475
2020 7.766635
2021 7.793272
2022 7.358050
ggplot(standard_dev_all_632_ad, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

Morrisey-corrected standard deviation of SNOTEL 632 average temperatures for water years 1986-2021

MK & SS 632 (corrected)

sd_mk_632_ad <- mk.test(standard_dev_all_632_ad$sd_2)
print(sd_mk_632_ad)
## 
##  Mann-Kendall trend test
## 
## data:  standard_dev_all_632_ad$sd_2
## z = 0.9518, n = 31, p-value = 0.3412
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##   57.0000000 3461.6666667    0.1225806
sd_sens_632_ad <- sens.slope(standard_dev_all_632_ad$sd_2)
print(sd_sens_632_ad)
## 
##  Sen's slope
## 
## data:  standard_dev_all_632_ad$sd_2
## z = 0.9518, n = 31, p-value = 0.3412
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01086011  0.02505935
## sample estimates:
## Sen's slope 
## 0.008387789